Phase diagrams and Thomas-Fermi estimates for spin-orbit coupled Bose-Einstein 

Condensates under rotation 



o 

u 

< 

00 



CZ2 

i 

c3 



i 

c 
o 
o 



> 

00 

o 



X 

S-H 



Amandine Aftalion 1 & Peter Mason 2 

1 CNRS & Universite Versailles-Saint-Quentin-en-Yvelines, 

Laboratoire de Mathematiques de Versailles, CNRS UMR 8100, 

45 avenue des Etats-Unis, 78035 Versailles Cedex, France 

2 Joint Quantum Centre (JQC) Durham-Newcastle, Department of Physics, 

Durham University, Durham DH1 3LE, United Kingdom. 

(Dated: April 9, 2013) 

We provide complete phase diagrams describing the ground state of a trapped spinor BEC under 
the combined effects of rotation and a Rashba spin orbit coupling. The interplay between the 
different parameters (magnitude of rotation, strength of the spin orbit coupling and interaction) 
leads to a rich ground state physics that we classify. We explain some features analytically in the 
Thomas Fermi approximation, writing the problem in terms of the total density, total phase and 
spin. In some regions of the phase diagrams, we relate the patterns to a ferromagnetic energy. 
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I. INTRODUCTION 

Bose Einstein condensates (BEC's) provide a unique 
experimental and theoretical testing ground for many 
macroscopic quantum phenomena. One such area that 
has recently attracted a lot of attention is the engineer- 
ing of synthetic non-Abelian gauge potentials coupled to 
neutral atoms [TH1] to create a spin-orbit coupled Bose- 
Einstein condensate [SJ [5] , where the internal spin states 
and the orbital momentum of the atoms are coupled. 
These spin-orbit coupled condensates support a variety 
of ground state density profiles; for instance in the most 
straightforward case of a spin-1/2 condensate [THTrJ]. the 
density either displays a plane wave or a striped wave. 
The transition between the two depends on the interac- 
tion parameters. However, if one is to consider a spin-1 or 
spin-2 condensate then more exotic ground state profiles, 
based on the helical modulation of the order parameter, 
can be created [T7H2U] . 

In addition to the various ground state density pro- 
files, one can look to the basic elementary excitations 
created in these spin-orbit condensates, like the vortex 
[2lTf2"3] . dark soliton [24], or bright soliton [25] in spin- 
1/2 condensates, or the skyrmion in spin-1 and spin-2 
condensates [26 29 . In contrast to a single component 
or two-component condensate, where the appearance and 
energetic stability of the elementary excitations are de- 
pendent on a rotation of the system to impart angular 
momentum, spin-orbit coupled condensates naturally im- 
part momentum through the coupling of the internal spin 
and orbital momentum of the atoms. But when combin- 
ing both the spin-orbit coupling and the rotation, vari- 
ous novel features have been predicted to occur [29TI33J . 
Through a suitable control of the condensate, an exper- 
imental scheme for rotating spin orbit coupled conden- 
sates has been proposed in [3D]. 

The aim of this paper is to study the combined effect 



of a Rashba spin-orbit coupling and a rotation on spinor 
BEC's for spin-1/2 condensates. The interplay between 
trap energy, spin-orbit coupling and interaction leads to 
a rich ground state physics: stripe phases, half vortices 
or vortex lattices with some behaviours reminiscent of 
vortex lattices appearing for fast rotating condensates 
[34l 135] . We provide a complete phase diagram accord- 
ing to the magnitude of rotation, spin orbit coupling and 
interaction. Some features have been analyzed by Sub- 
hasis et al. [23], Zhou et al. [32] and Xu k. Han [35] . 
but here we want to investigate a full phase diagram be- 
haviour. 

Our paper is organised as follows. In Section II we in- 
troduce the energy functional in terms of individual wave 
functions before making the transformation to the non- 
linear Sigma model where the energy is instead written 
in terms of the total density and a spin density. In Sec- 
tion III we provide numerically determined phase space 
diagrams for the ground states of the condensate as func- 
tions of the rotation, spin-orbit coupling and interac- 
tion. We explain some features analytically by using a 
Thomas-Fermi approximation in Section IV. 



II. PROBLEM STATEMENT AND ENERGY 
FUNCTIONAL 



We are interested in a two-dimensional (x, y) rotating 
spin-coupled Bose-Einstein condensate. This has the fol- 
lowing non-dimensional energy functional in terms of the 



wave functions ipx and %p2 
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under the constraint that JIV'il 2 + \ifa | 2 = N, N be- 
ing the number of atoms. Here, gk is the self interac- 
tion of each component (intracomponent coupling) that 
we will later take to be equal for both components, 312 
measures the effect of interaction between the two com- 
ponents (intercomponent coupling) and Q is the rota- 
tional velocity, applied equally to both components, with 
L z = —i(xd y — yd x ) the angular momentum operator act- 
ing in the z direction. We consider a Rashba spin-orbit 
interaction strength, k being equal in both the x and y 
direction. 



One of the key ingredients in the analysis will be to use 
the nonlinear Sigma model introduced for two compo- 
nent condensates in the absence of a spin-orbit coupling 
[36H38] ; that is to write the energy in terms of the total 
density p, 



IV'il' + l^l 5 
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and a normalised complex- valued spinor x — [xitX2] T - 
the wave functions can be decomposed as ipi = y/pxit 
■02 = \[PX2 where |xi| 2 + IX2I 2 = 1- We define the spin 
density S = X&X, where <x = {<Jx,&y,<Jz) are the Pauli 
matrices, with the components of S following as 



Sx =XlX2+X2Xl> 


(3a) 
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(3b) 
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such that \S\ 2 = 1 everywhere. For a rotating conden- 
sate, it is natural to introduce 
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where O — 81 + 02, Ofc is the phase of ipk, that is ipk = 
yfp\Xk\e mk and R = S y VS x - S x VS y . This allows us to 



rewrite the energy functional (TTT) as 
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A derivation to this form of the energy from (nl) and to 
other forms is given in the Appendix. 

Note that our choice of the effective spin-orbit Hamil- 
tonian [Eq. (fl])] is assumed to remain stationary in the 
rotating frame. This is in contrast to the experimental 
schemes proposed in [30] . where there is a time depen- 
dence inherent in the Hamiltonian. We justify our as- 
sumption and use of a time-independent Hamiltonian on 
two fronts: firstly the probable small effect that the time 
dependent terms will have on the ground state (we note 
in particular Fig. 1(c) of Ref. [30 in which a regular 
vortex lattice is present in both components for a large 
spin-coupling and a relatively large rotation); secondly, 
the need to perform meaningful analytical analysis on the 
ground state profiles requires a 'from principles' approach 
whereby only the fundamental terms of the Hamiltonian 
are considered, that is the spin-orbit coupling, the rota- 
tion and the interaction, as written in the Hamiltonian 
u\. Furthermore to this last point, the experimental in- 
frastructure to create a rotating spin-orbit condensate is 
relatively new, and there remains the possibility that a 
new experimental scheme that fully justifies the use of 
a time-independent Hamiltonian could be proposed. To 
this end, we believe that our phase diagrams provide in- 
teresting and relevant information on the ground states 
of the rotating spin-orbit coupled condensates. 



III. DESCRIPTION OF THE PHASE DIAGRAMS 

We wish to describe the ground state wave functions 
of the spin-orbit condensate. In what follows, we assume 
5i = 92 = 9 and set S = 9vil 9-, which measures the effect 
of interaction between the two components. The exper- 
iments of [HI IS] have gN large (the Thomas- Fermi limit). 
Therefore, our analysis will also be in the case gN large 
and our system is then described by three parameters: 



f2, the rotational velocity; k, the spin-orbit interaction 
strength; and 8. 

We first present numerically obtained phase diagrams 
for these three parameters, with a Thomas-Fermi analysis 
following in the next section. These simulations are con- 
ducted on the coupled Gross-Pitaevskii equations that 
result from the energy functional (nj) through the varia- 
tion idijjk/dt = 5E/6i[)* k for fc = 1, 2: 
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We simulate in imaginary time using the following values 
of parameters: g — 4 and N — 200, together with f2 £ 
[0,1), k > and 8 > 0. These parameters place us 
in the Thomas-Fermi regime. For each parameter set we 
classify the ground state according to the densities, \ipk\ 2 , 
and the spin densities, S. In general, it is difficult to find 
the true minimizing energy state. But the use of various 
initial data converging to the same (or similar) final state 
allows us to determine that the true ground state will be 
of the same pattern as the one that we exhibit. We break 
our analysis into three sections; Q — 0, fi small and f2 
large. 
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FIG. 1: Density profiles (left column, component-1 and right 
column, component-2) when Q. = and (5, k) equals (I) 
(0.25,1), (II) (1.44,0.5) and (III) (1.44,2). 



A. n = o 

We begin by considering the non-rotating spin-orbit 
condensate in which the active parameters are k and 8. 
In the case when n — 0, we are left with a two-component 
condensate coupled exclusively by the intercomponent in- 
teraction strength related by 8. In this case, the ground 
state density profiles of the condensates are either, for 
8 < 1, two co-existing disks (of equal radii), or for 5 > 1, 
one of the components is a disk while the other is identi- 
cally zero [37]. For all 8, there are never any topological 
defects created in the condensate. 

Turning on the spin-orbit coupling term so that k^0 
provides a system which has recently been considered 
in the literature by a number of authors [5HTUI [T71 IH)I 
[2"2l |2"3"] . When 8 < 1, the two-components remain co- 
existing and disk-shaped for all k (Fig.JTVl)). If, on the 
other hand, 8 > 1, for small k (Fig. flfll)), then one 
component is a disk, in which most of the particles reside 
and is surrounded by a thin, low populated annulus for 
the other component. Note that which component is a 
disk and which is an annulus cannot be determined. As 
k gets larger, the two components are disks with stripes 
present in the density profiles (Fig. ITTlII) ) : these density 
profiles were studied in [51 [53] . We include a summary 
figure detailing these three ground state profiles in Fig. 



m 

When 5 < 1, we never see any topological defects in the 
density profiles of the coexisting disk shaped condensates. 
However, singularities can be seen in the total phase pro- 
files 0, where 8 = 0i -I- 02 and the ©& are the phases 
of the ipk (the individual phase profiles, ©i and 02 also 
contain singularities). For instance in Fig. [2] we plot 
0, dQ/dx and d<d/dy for two different parameter sets, 
(8, k) = (0.25,0.5) and (0.25,1). For k = 0.5, there are 
no singularities present within the Thomas-Fermi radius 
(see Fig. [2la,I)), however for n — 1 we see the develop- 
ment of singularities (see Fig. [2jb,I)). Notice that (away 
from the singularities) dQ/dx is a (negative) non-zero 
constant (Fig. [^II)) whereas dQ/dy = (Fig. [2J1II)) 
when n takes either of these values. A further increase in 
k introduces more singularities into the total phase; see 
Fig. [3] where we plot the profiles of dQ/dx and dQ/dy 
for (8,k) = (0.25,4.75). In fact, here (d/dx + d/dy)Q/2 
and (d/dx — d/dy)Q/2 play the same role as dQ/dx and 
dQ/dy do for smaller k. 

In addition to the development of singularities in 0, 
we must look to the profile of S. The spin density vec- 
tor exhibits two different states: the first being the case 
where S x « 1 and S y « S z w everywhere, except for a 
thin layer near the Thomas-Fermi surface (see Fig. Ba) 
where (8, k) = (0.25,0.5)) and the second being the case 
where S x ~ ±S y ~ l/v2 and S z « everywhere, again 
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FIG. 2: Phase plots (frame (I) 6, (II) dQ/dx and (III) dQ/dy) 
for numerical simulations carried out when Q. — for (a) 
(8,k) = (0.25,0.5) (left column) and (b) (S,k) = (0.25,1) 
(right column). Notice that (away from singularities) dQ/dx 
is a non-zero constant whereas dQ/dy is zero. The Thomas- 
Fermi radius is plotted (dark black circle), calculated in Eq. 
(12a). 





FIG. 3: Phase plots (frame (I) dQ/dx and (II) dQ/dy) for 
numerical simulations carried out when Q. = for (S, k) = 
(0.25,4.75). The Thomas-Fermi radius is plotted (dark black 
circle), calculated in Eq. (12a I. 



except for a thin layer near the Thomas-Fermi surface 
(see Fig. Em) where (S, k) = (0.25,4.75) and in this par- 
ticular example S x « +S y ). 

The transition between these two states occurs for 
K ~ 3.5 (with S — 0.25) and is independent of the de- 
velopment of the first singularity in the total phase pro- 
file. Thus the cases considered in Fig. [2] both have 
S x w 1 and S y « S z ~ 0, whereas the case in Fig. [3] 
has S x « ±S y « l/\/2 and S z w 0. We will relate this 



FIG. 4: Spin density plots (frame (I) S x , (II) S y and (III) 
S y ) for numerical simulations carried out when fi = for (a) 
(5,k) = (0.25,0.5) (left column) and (b) (6,k) = (0.25,4.75) 
(right column). The Thomas- Fermi radius is plotted (dark 
black circle), calculated in Eq. ( 12a I. 



behaviour to ferromagnetic problems in Sect. IV. 



B. fi small 

We now proceed to the case where $7 and k are in gen- 
eral non-zero. We first present a k — S phase diagram for 
fl = 0.1 (small rotation) in Fig. [5] in which four distinct 
regions are present. We identify these as (i) two disks 
with no defects, (ii) two disks with domains, (iii) a single 
component and (iv) stripes. Each region on the phase di- 
agram of Fig. [5] contains a sample density profile from a 
simulation carried out within that region (the simulation 
parameters are noted in the figure). The key difference 
between a phase diagram with Q = and f2 small is 
the development of region (ii), which is not present when 

n = o. 

Along the k — axis, we revert to the case of a rotating 
two-component condensate for which the ground state 
profiles are two co-existing disks (S < 1) or there is spatial 
separation of the components - one component is a disk 
and the other has a zero wave function (6 > 1) [37|,l38j. 
These behaviours are still present for k small (regions (i) 
and (iii) respectively of Fig. [5| . The profiles of the spin 
densities related to these two regions are straightforward: 
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FIG. 5: k — S phase diagram with Q — 0.1. The numer- 
ical parameters are taken as g — 4 and TV = 200. There 
are four identified regions: (i) two disks with no defects, (ii) 
two disks with domains, (iii) a single component and (iv) 
stripes. Each region has a typical density plot for each compo- 
nent (left panels, component-1 and right panels, component- 
2). The numerical values of these simulations correspond to: 
(i) (5,k) = (0.25,0.5); (ii) (0.5,5); (iii) (1.69,0.25) and (iv) 
(1.44,4). We analyse regions (i) and (ii) in more detail in 
Fig.'sp]and[7] The transition between regions (iii) and (iv) 
is described in Sect. HID. 



in region (i) we have S = (S x , S v , S z ) « (1,0,0) - much 
the same as in Fig. Rla), whereas in region (iii) we have 
S « (0,0, 1). As k becomes larger, modulations of the 
density profiles begin to occur. For 8 > 1, the spatial 
separation of the two-components is lost. This occurs 
for relatively small values of k (see Fig. [51 where the 
boundary between regions (iii) and (iv) occurs for < 
K < 1 for 1 < S < 2, and Sect. HID). Over this boundary, 
a stripe pattern emerges in both components. 

For S < 1, the co-existing disk-shaped components 
each develop vortices that arrange themselves along 
bands of each component. We classify this region as 
the region in which both components are 'two disks with 
domains' (region (ii) in Fig. pj. The domain we refer 
to here is related to the profile of S. We notice - see 
Fig. [6] - that across some bands, the behaviours of the 
S x and S y components of the spin density change sign. 
For example, Fig. |6jlII,IV) plots the S x component and 
the S y component for the parameters (5, n) — (0.5, 1.25) 
and (S, k) — (0.5,5) [O = 0.1]. In this simulation with 
(S, k) — (0.5, 1.25), two bands of vortices have been cre- 
ated along the x axis, while for (5, k) = (0.5, 5) there are 
four bands of vortices, each along one of the principal 
axes. In both cases we see that for y > (< 0), S x > 
(< 0) and for x > (< 0), S y < (> 0). This creates do- 
mains within the S x and S y component profiles (we note 
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FIG. 6: Density plots (frame (I), component-1, and (II), 
component-2) and spin density plots (frame (III), S x , frame 
(IV), S y and frame (V), S z ) for region (ii) of the k — S with 
fi = 0.1 phase diagram of Fig. [5] The numerical simulations 
are carried out for (a) (5, k) — (0.5,1.25) (left column) and 
(b) (S,k) = (0.5,5) (right column). 



that S z ~ away from the vortex lines). For this par- 
ticular example, we say that there are two domains. A 
particular feature of the domain structure of the S x and 
S y is that, away from the vortex lines, they become (ap- 
proximately) constant. For example, in Fig. [6FlII,IV), 
S x w l%/2 (w -lV2) for y > (< 0) and S y w -ly/2 
(w l\/2) for x > (< 0) [note that Si + S^ w 1 as we 
have S z sa everywhere]. As k is increased to higher 
values, we see examples with more domains. 

When k = 0, the critical value 6 = 1 divides the con- 



densate profiles between co-existence and spatial sepa- 
ration. For the small value of 17 chosen for the phase 
diagram in Fig. [5j 5 — 1 as a critical value is still rele- 
vant. However as k increases, our numerical simulations 
indicate a blurring of the boundary between regions (ii) 
and (iv). To indicate the uncertainty in the location of 
this boundary for high k, we have used a dashed line in 
Fig. [5] above some arbitrary k. 

Another important feature of the ground states is 
the profile of the total phase, 9. This constitutes, to- 
gether with p and two components of S [i.e. two from 
(S x , S V ,S Z )], a complete reformulation (four variables) of 
the original problem in terms of the two wave functions. 
The total phase, and its gradient, are plotted in Fig. [Tjfor 
two numerical values, (S, ft) = (0.25,0.5) (corresponding 
to region (i)) and (S, ft) = (0.5,1.5) (corresponding to re- 
gion (ii)) with $7 = 0.1. When we are in region (i), where 
the density profiles display two disks with no vortices, the 
total phase (within the TF radius) varies linearly with x 
and shows no variation in y, as can be seen by frames 
(a,I-III), where in (a,II) we plot dQ/dx and in (a, III) we 
plot dQ/dy. The value of dQ/dx is constant and equal to 
— 1.01, whereas the value of dQ/dy is zero. At this stage 
we merely note that —1.01 ~ —2kS x , and ~ —2kS v , re- 
lationships that we will formalise later. In contrast, when 
we are in region (ii), where the density profiles display 
two disks with bands of vortices, the total phase displays 
singularities within the TF radius (the singularities orig- 
inate at the location of the vortices in each component); 
see Fig. [7fb,I). Of more interest are the plots of dQ/dx 
(b,II) and dQ/dy (b,III). In both these cases, away from 
the singularities, the x- or y-derivative of the total phase 
is approximately constant, and odd, respectively, with 
the x and y axis. For example, dQ/dx ~ —3 for y > 
and dQ/dx ps 3 for y < 0. On the other hand, dQ/dy is 
small everywhere but the singularity lines are odd with 
respect to the y axis. Again this is in agreement with the 
relationship (dQ/dx, dQ/dy) — — 2n(S x , S y ) at leading 
order. 



C. Large fl 

If instead we look to the large $7 limit, then we see the 
rotational effect dominating. Figure [8] shows a k— 6 phase 
diagram for 17 = 0.9 (large rotation: note that 17 must 
stay below max(17) = 1) in which three distinct regions 
are present. We identify these as (iii) a single component 
(as in Fig. [5]) , (v) two disks with vortex lattices and (vi) 
two annuli with vortex lattices. Each region on the phase 
diagram of Fig. [8] contains a sample density profile from a 
simulation carried out within that region (the simulation 
parameters are noted in the figure). 

Again if we consider the k = axis then we revert 
to the two-component condensate rotating at high an- 
gular velocities [38 . In these cases, the large rotational 
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FIG. 7: Phase plots (frame (I) 6, (II) dQ/dx and (III) dQ/dy) 
for numerical simulations carried out when Cl — 0.1 for (a) 
(5,k) = (0.25,0.5) (left column) and (b) (5,k) = (0.5,1.5) 
(right column). Notice that (a,II) is a non-zero constant 
whereas (a, III) is zero and that (away from singularities) 
(b,II) is a non-zero positive (negative) constant for negative 
(positive) y whereas (b,III) is zero. The Thomas-Fermi radius 
is plotted (dark black circle). 



effect leads to angular momentum being imparted onto 
the condensate and therefore to the existence of vor- 
tices. Thus the condensate is either comprised of two 
co-existing disk-shaped components both with a triangu- 
lar (coreless) vortex lattice (for 8 < 1), or of a single com- 
ponent with a triangular vortex lattice (for S > 1, there 
is spatial separation where the other component has zero 
wave function). As the spin-orbit-coupling interaction 
is included into the problem (ft becomes non-zero), then 
the spatial separation region ((iii) in the phase diagram of 
Fig. [8]) is soon lost. Thus, for high rotation, the inclusion 
of a (small) spin-orbit-coupling interaction strength leads 
to the spatially separated components becoming equally 
populated (i.e. there is a transition from S « (0, 0, 1) to 
S = (S x , S y , 0) where the S x and S y components of the 
spin density are in general non-zero) . Figure K)ta) plots 
the component densities and spin densities for region (v) 
of the 17 = 0.9 phase diagram. 

It is interesting to note the difference in the spin den- 
sity profiles of this high rotation, low spin-orbit coupling 
case to the low rotation, high spin-orbit coupling case of 
Fig. [6| particularly the form of S x and S y . While in 
Fig. |6]the S x and S y are almost constant, in Fig. [9] the 



1.4 
1.2 

1 
0.8 
0.6 
0.4 
0.2 

0, 



2 annuli with vortex lattices 






2 disks with vortex lattices 







Single component 



0.5 




FIG. 8: k — 8 phase diagram with Q — 0.9. The numer- 
ical parameters are taken as g — 4 and N = 200. There 
are three identified regions: (iii) a single component, (v) two 
disks with vortex lattices and (vi) two annuli with vortex lat- 
tices. Each region has a typical density plot for each compo- 
nent (left panels, component-1 and right panels, component- 
2). The numerical values of these simulations correspond to: 
(iii) (5,k) = (1.25,0.025); (v) (0.25,0.5) and (vi) (0.25,1.25). 
The boundary between regions (v) and (vi) is plotted accord- 
ing to the numerical simulations (solid line) and analytically 
(dashed line), calculated according to Eq. ( 34[). We analyse 
regions (v) and (vi) in more detail in Fig. |9 



S x and S y are sine/cosine-like functions. We will show 
this to be the case later, but we note for now that, in 
essence, we see a smooth sine/cosine-like function for S x 
and S y in the rotation dominating regime, whereas in the 
spin-orbit dominating regime we see S x and S y becom- 
ing constants with sharp transitions over boundary lines 
(that correspond to the lines of vortices and the defini- 
tion used in this paper for the domain boundary). The 
vortices of each component correspond to singularities in 
the S x , S y , S z components: pairs of upwards and down- 
wards spikes. 

The importance of the 6 = 1 value is soon relegated 
as re becomes non-zero. Instead, for all S and k > 0.05, 
the components are equally populated and display, in all 
cases, triangular vortex lattices. However the geometries 
of the components can change, dependent upon re and S 
(and Q), although both components will, for any given 
parameters, have the same geometry In particular, one 
can see the development of two annular components as re 
is increased. These annular components still preserve the 
vortex lattices, and the sine/cosine-like form of the spin 
density - see Fig. [9lb). In the next section we find an 
analytical expression for the critical parameters at which 
the geometry changes from two disks to two annuli. This 
analytical result (dashed line) can be compared to the 




FIG. 9: Density plots (frame (I), component-1, and (II), 
component-2) and spin density plots (frame (III), S x , frame 
(IV), S y and frame (V), S z ) for regions (v) and (vi) of the 
K — 8 with Q, — 0.9 phase diagram of Fig. [8] The numeri- 
cal simulations are carried out for (a) (8, k) — (0.25, 0.5) (left 
column) and (b) (0.25, 1.25) (right column) respectively. 



numerical simulations (solid line) in Fig. ^j 



D. Transition from spatially separated components 
to stripe phases 

In Fig. [l|lI,III) we show the density profiles for a non- 
rotating condensate with S > 1. When re is small the 
density profiles are a disk plus annulus (TTlI)] and when 
re is 'large' the density profiles are stripes fijll!)]. If we 










FIG. 10: Density plots (left column, component-1 and right 
column, component-2) with £1 = 0.1 and 8 — 2. The numer- 
ical simulations are carried out for (I) k = 0, (II) k = 0.25, 
(III) k = 1.25 and (IV) K = 1.5. 



FIG. 11: Density plots (left column, component-1 and right 
column, component-2) with £1 — 0.5 and 5 = 2. The numer- 
ical simulations are carried out for (I) k = 0, (II) k = 0.25, 
(III) k = 0.75 and (IV) k = 1.75. 



look solely at S > 1, then the same features are present 
when O = 0.1. In the phase diagram of Fig. [5j we 
have drawn the transition between the two regimes [from 
regime (iii) to regime (iv)] as being instantaneous, but 
in reality there is a smooth transition from one profile to 
the other. In Fig. [lO] we show some density profiles that 
correspond to values of k taken around this transition. 
For larger rotation the main predominant transition is 
from the spatially separated components to a two disk 
state. In Fig. [II] we also show some density profiles that 
correspond to this transition. 



E. Q-5 phase diagrams 

We have up to now only presented k-S phase diagrams 
with the value of the rotation held constant (either fl — 
0.1 for Fig. [5|or Q = 0.9 for Fig. [8]). In addition to these 
we present two fl-5 phase diagrams with the value of k 



held constant: in Fig. [12(a) we take k = 1 (small) and 
in Fig. pFb) k = 8 (large). 



IV. ANALYSIS IN THE THOMAS-FERMI 
REGIME 



In this Section we perform an analysis based on the 
Thomas-Fermi approximation to describe various fea- 
tures of the phase diagrams that we introduced in the 
previous section. In particular, we will concentrate on the 
case S < 1, that is on the symmetry preserving ground 
states (those featured in Fig. [IJa) and regions (i)-(iii), 
(v) and (vi) of the phase diagrams shown in Fig. 's [5] and 
[8]). Our starting point is the energy functional [Eq. pi] 
written in terms of the total density and the spin density, 
for which we consider, as in the numerical simulations, 
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fi 



fi 
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2 annuli 
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Stripes 



2 disks no defects 
(15 



1.5 



terms in Eq. &. We divide our analysis into looking at 
the cases of zero rotation and non-zero rotation which we 
further divide into low and high rotation. 



A. No Rotation 

Let us begin by looking at the case where there is no 
rotation, £1 — 0. Since 8 < 1, that is C2 > 0, we can 
assume in the energy (|8j), that C2S 2 is negligible in front 
of Cq, and v c r ~ V0/2 which leads to 



E = 



/i(V^f 



lO> 2 



kS-VxS 



^(vef + ^-ve + V + co^ dV 



(9) 



This energy leads to two orders of magnitude, one for p 
and the other for S and O. Thus, we have to separately 
minimize 



and 



E P = 



Es,e = 



\ (vVp) 2 



P 2 

r 



c - 



p 



d z r, 



(10a) 



l\(vsf 



J(V6) 2 



iS-V x S 



KS ± Ve) d 2 r. (10b) 



FIG. 12: Q — 8 phase diagrams with (a) n = 1 and (b) k — 8. 
The numerical parameters are taken as (? = 4 and N = 200. 
There are six identified regions: (i) two disks with no defects, 
(ii) two disks with domains, (iii) a single component, (iv) 
stripes, (v) two disks with vortex lattices and (vi) two annuli 
with vortex lattices. The boundary between regions (v) and 
(vi) is plotted according to the numerical simulations (solid 
line) and analytically (dashed line), calculated according to 
Eq. pill. 



gi = 92 which gives c% = 0: 

E=j\ {VVP? + § (VS) 2 + P - (v cS - fl x r) 



pn[ S±-v cS + -S- Vx S 



+ ^l-tfy + {c a + c 2 Sir- cfr. 



(8) 



Under the assumption that gN is large, we are in the 
Thomas-Fermi limit, which allows us to make various ap- 
proximations regarding the importance of the individual 



Eq. (10a) yields the Thomas-Fermi profile for p that 



we will discuss below. Eq. (10b) leads to two coupled 
problems for 8 and S. 

Thomas-Fermi Density Profiles. In the Thomas-Fermi 



approximation, the minimization of (10a) yields 



P 



1 
Co 



/' 



-r 

2 



(11) 



where p is the chemical potential. This Thomas-Fermi 
density profile has a harmonic trapping potential and so 
the components will always be disk shaped. To complete 
the analysis we use the normalisation condition to obtain 



R = 



P = 



/4JVc \ 1/4 

(^ 



(12a) 
(12b) 



where R is the Thomas-Fermi radius. This value fits very 
well with the numerical computations and has been used 
to draw the limiting region of the condensate in Fig.'s 

da 
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Equations for O and S. The minimization of ( 10b ) 
leads to an equation for written as 



B. Non-zero Rotation 



V ■(p{VQ + 2nS±)) =0. 



(13) 



1. Low Rotation 



This is reminiscent of the continuity equation written in 
|39| . For small k, S± can be written as a gradient so 
that 



V6 + 2k5_l = 0. 



(14) 



The numerics indeed illustrate that S x = 1, S y = S z = 
and that V0 = —2kS± everywhere (see Fig.'s [2Fa) 
andWa)). Therefore, when E g. ([l4| is satisfied, the 
minimization of Es,&, given by (10b), becomes 



E 



s.» = I f (\ (VS) 2 + K S-X7xS + n 2 (S 2 - l)\ d 2 r, 

(15) 
since Sj_ = 1 — S 2 . The ground state of this energy for 
small k should be proved to be close to S = (1, 0, 0). 



When the rotation is included into the problem then 
the energy ^ also leads to decoupled problems, one in 
p, and one for S and O: 



Er = / \ Wf>? 



P 2 

-r 

2 



c - 



P 



d 2 r, 



(19a) 



E. 



s,e = 



jd\(WS) 2 + K S-WxS + K 2 S 2 

+ J (V6 + 2kS±) 2 -Ve-fixr (S§h) 



The Thomas-Fermi expression for p does not change with 



respect to the O = expression [Eq. (11 1]. We still have 



that (14 1 holds except on singularity lines. New singu- 



larities related to the rotation emerge, on the boundary 
of the domain regions of S± as illustrated in Fig. [6j 



As k is increased, then Eq. (14 1 cannot be preserved: 



it is satisfied everywhere except along singularity lines: 
see Fig. [2lb). Indeed, V x S± is no longer zero but 
approaches a measure supported on the boundary of the 
Thomas-Fermi disk, which produces singularities in Eq. 
pT] ). If k is further increased, then the behaviour of 
S changes and S - (l/%/2, 1/^2, 0) as in Fig. gb). 
Fig [3] plots the phase profiles associated with this case 
and we see in this example, away from the singulari- 
ties, we have {dQ/dx, dQ/dy) = (-6.674,-6.674), that 



is — 2kS±. Eq. (13) implies the existence of a function £ 
such that 



p(V9 + 2kS , _ l ) = V-^, 
where V^ = (—d y ,d x ). This implies 



(16) 



2. High Rotation 



When the rotation is increased, the energy becomes 
(see Appendix) 

E = l\ ( V v^) 2 + f { y S f + \ ( v -« + kS ± -Ux r) 2 

1 



pQ,n{—yS x + xSy) + -puS -VxS 



+ l -pn 2 S 2 z + f (1 - n 2 )r 2 + (co + c 2 S 2 z ) P - d 2 r. 

(20) 



V- [~V£ 



2«;VxS 1 -V-(V 1 e). (17) 



If V x S± is a measure supported on the boundary of 
the Thomas-Fermi disk, then V0 should have singularity 
lines. This leads to a minimization problem for S which 
is related to ferromagnetic equations for thin films: 



1 

2 V4 



2o2 



(VS) + kS ■ V x S + k 2 S 



1 
8p 



|V£| 2 d\ 

(18) 
where £ is the solution of (16). For intermediate k, this 



leads to boundary singularities for S± and S z 



In the Thomas Fermi regime, this decouples into a prob- 
lem for p, (for which the kinetic energy is negligible and 
C2S 2 can be neglected in front of Co) and a problem for 



v e f[ and S. This allows us to rewrite Eq. (20) as follows 



E n = f P QK(-yS x + xS y ) + |(1 - f> 2 )r 2 + c 



P d 2 r, 
2 

(21) 



E vm ,s = J I (VS) 2 + I (v cS + K S ± -£lx rf 

+ -pnS ■ V x S + -pn 2 S 2 d 2 r. 
2 2 



(22) 
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The minimization of (21) leads to two Euler-Lagrange 
equations: 

-(1 -VL 2 )r 2 + Ktt(-yS x + xS y ) + c p = n, (23a) 

K %Ti(l-^r 1/2 4) = 0. (23b) 

We note that, at leading order, S 2 + S 2 ss 1 , so that Eq. 
(23b) gives 



S~ (sin 9, -cos 6»,0). (24) 

The next order in S has to be analyzed with the mini- 



mization of (22), which will lead to the vortex lattice. 



Substituting (24) into Eq. (23a) gives 

p= — ( fi--(l- n 2 )r 2 + nVtr 
co V 2 



(25) 



We can use the normalisation condition to find /i. How- 
ever first we must consider the two possible geometries: 
by our assumption that S z « away from the defects, we 
expect the two components to share the same geometry. 
The numerical simulations presented in Fig. M (and also 
present in [33 ) indicate that both components are either 
disks or annuli. 

Two disks. When the components are both disks (we 
show later that this corresponds to /i > 0), then we can 
use the normalisation condition, integrating from r = 
to the outer boundary, r = R 7 to find that R is given as 
the solution of the quartic 



I? 4 



Ann 



-R J 



ANc 



3(1 -fi 2 ) 7r(l-fi 2 ) 

which then gives the chemical potential as 

/!= -(l-fi 2 )i? 2 -KflR. 



0. 



(26) 



(27) 



Two annuli. When the components are both annuli, 
we can analyse an effective potential, defined from Eq. 
([25) as V e (r) = |(1 - fi 2 )r 2 - nVtr. The total density is 



zero when this effective potential is equal to the chemical 
potential: i.e. V e (r) = \x. Thus 



1 



/'■ 



(1 -Q 2 )r 2 -KQr 



K n± J n 2 n 2 + 2(1- n 2 )fi 

=> r = 

(i-fi 2 ) 

from which we identify 

an- y/K 2 n 2 + 2(i-n 2 )n 



Ri = 



R2 



(i-n 2 ) 



K,n + y/ K 2 n 2 + 2(i-n 2 )[i 
(i-fi 2 ) ' 



(28) 

(29a) 
(29b) 



with the inner radius i?i of the annulus and the outer 
radius i? 2 - For i?i to exist we must have \x < 0, i.e. 
fi > implies that both components are disks and /1 < 
implies that both components are annuli. We thus here 
assume that /i < 0. Note that 



Ri + R% — 
i?2 — R\ = 



2kQ 

(1-ft 2 )' 

2y / K 2 n 2 + 2(i-n 2 )n 



(1-n 2 ) 

We can find /i from the normalisation condition as 



(30a) 
(30b) 



1 



M = 



-K 2 n 2 i mc o0--sr) 



2^3\2/3^ 



2(1 -n 2 ) v V 47rKfi 

which gives R\ and R2 explicitly as 



(31) 



k£1 



yK c \l c ) 



Rl - (T^) " [nnii'-k) 3 ) (32a) 



R 2 = 



kQ 



yK c \L c ) 



, (32b) 



(1-n 2 ) V K ^(!-^ 2 ) 3 , 

which leaves the width of the annulus, d = R2 — Ri, as 

(33a) 



where 



a = 2 



(KcOc) 4/3 



(33b) 



is a constant, with the k c and Q c critical values of k and 
n (and 512) that can be found from setting (i = and 
using the normalisation condition. This gives 



yK c \l c ) 



3Nc 

47T 



(34) 



Comparison to Phase Diagrams. A phase diagram iden- 
tifying the numerically determined ground states for the 
case of large rotation is given in Fig. [S] We see the ap- 
pearance of three regimes, of which the above analysis 
pertains to the regimes in which there exists two disks 
with vortex lattices and two annuli with vortex lattices. 
The boundary between these two regimes has been cal- 
culated analytically to be given by Eq. (34), which we 



include on the phase diagram of Fig. [8] (dashed line), to- 
gether with the numerically determined boundary (solid 
line). One can see a good agreement between the theory 
and numerics. 



The width of the annulus given by ( 30b ) becomes thin 



ner as the product k£1 increases. This scenario is remi- 
niscent of the rotating single component condensate that 
is trapped by a harmonic plus quartic trapping potential 
[40] in which an annulus develops as the rotation is in- 
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creased, and in which the width of the annulus becomes 
smaller. In the case of a condensate held by a harmonic 
plus quartic trapping potential, there is no upper limit 
to the rotation since the quartic term acts to keep the 
condensate bounded for all Q. This increasing rotation 
leads to the development of a giant vortex (a large circu- 
lation) inside the annulus. At the same time the width 
of the annulus is decreasing such that the condensate can 
no longer support any vortices in the condensate bulk. A 
similar situation can occur if one considers a harmonic 
plus Gaussian trapping potential (a toroidal trap), al- 
though in this case the width of the annulus is dependent 
on additional factors, notably the strength and 'waist' of 
the Gaussian term (which is generally taken to be cen- 
tred at the origin) and an upper limit on the rotation 
which must be enforced to ensure the condensate stays 
bounded [ITj . 

For the spin-orbit condensate that we consider in this 
paper, the width of the annulus is determined by both the 
rotation frequency, fi, and the spin-orbit coupling term, 
k. While the rotation frequency is again bounded (its 
upper limit is 1), there is no limit to the value that k can 
take (at least theoretically). As such, for a sufficiently 
large k, a thin annulus with a large circulation can be 
created (as has been seen in Fig. 2 of Ref. [33]). 



Appendix: Derivation of the Energy Functional (pil) 



1. First formulation 



We show in this Appendix how to derive the energy 
functional 



E = 



\Wp) 2 + ^(vs) 2 + f(« off -n x rf 
1 



pn[S ± -v cff +-S-V x S 



+ '-(i-n z y + (c Q + c 1 s z + c 2 si) 



%\P 



d 2 r, 



(35) 



in terms of the non-linear Sigma model. We start with 
the energy functional given in terms of the wave functions 
ipk [Eq. IT], rewritten here 



E = 



/ E [l\Wk\ 2 + ~r i \il>k\ a -WkLzil>k 

J fc=l,2 \ 



+ fw 



K^fc 



+ 312 1 ^1 1 2 1^2 



\ dlp3-k + , 1 . 3 _ fc gV>3-fc 

dx dy 

d 2 r. 



(36) 



V. CONCLUSION 



We have provided phase diagrams in terms of the mag- 
nitude of rotation, the strength of the spin orbit coupling 
and interactions. We have found that plotting the total 
phase and the components of the spin leads to an in- 
teresting classification of the ground states. In the case 
of coexisting condensates, the Thomas Fermi approxima- 
tion for the total density leads to a simplification of the 
energy. We are able to determine the boundary between 
regions of disks and annuli leading to vortex lattices at 
high rotation, and to derive a ferromagnetic energy en- 
compassing the singularity behaviours for the total phase. 



In [38 , we showed how to transform the energy functional 
of the rotating two-component condensate (i.e the above 
energy functional with k = 0) into one given in terms of 
the total density p, the total phase O and the spin density 
S. So we split the above energy functional into terms 
independent of the spin-coupling and terms dependent 
on the spin-coupling, E = E t + E so , where [55] 



E,= 



E ui v ^i 2 + \am 2 - nr k L,i>k + f m' 

fc=l,2 ^ 
f3l2|-0l| 2 |^2| 2 



d 2 1 



\{V^-p) 2 + P -{VS) 2 + P -{v c ^nx.r) 2 



P 



+ ^i-n z y + (c a + Cl s z + c 2 si) 



2\P 



d 2 r, 



(37) 
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and 

E sn : 



J k E ^ 



fe=l,2 



.dip 3 -k 
dx 



-1) 



3_ k dlp a - k 



dy 



d 2 r. 



(38) 
We work with Eq. (38), and use the non-linear Sigma 
model where we have tpk = \fPXk, where p — |V'i| 2 + |' , /'2| 2 
and the \k are related to S by Eq.'s pp. Note that 
l^l 2 = 1 everywhere, so that one of the components of S 
is given in terms of the other two. 
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Upon substitution of ipk — y/pxk in to Eq. (38 1, we 
find that 



*~ = -/f(i-s2) 1/a E 



e i(63-fc-6fc) x 



fc=l,2 



dx 

(-i)V 



as, 

2(l + (-l) fc S 2 ) 2 V" dx 



+ (-1) 
rf 2 r 



a-* as. 



<9y 



(39) 



where we have written the spinor x& in terms of its am- 
plitude and phase: \k — |Xfe| ex P(*Qfe)- Furthermore, we 
use the identities S z = |xi| 2 -|X2| 2 and |xi| 2 + IX2| 2 = 1 
togive| Xfe | 2 = (l + (-l) 3 - fc ^)/2. 



Next we note that S x = 2|xi 1 1x2 1 cos(6i — 62) = (1 
S 2 ) 1 / 2 cos(0i - 6 2 ) and S, = -2| Xl || X2 | sin^x - 2 ) 
-(1 - S 2 ) 1 / 2 sin(9i - 6 2 ), which allows us to write 

-i{l-Sl) 1 ' 2 J2 e 4 ( e — - 0fc )x 

fc=l,2 



.ae a _ 



3-fc 



dx 



•1) 



3 _ k dQ 3 - k 



dy 



sl^(e 1 + e 2 ) + ^(e 1 -e2) 



s(^(e 1 + e2) 



d 



(9! - e 2 



Si -ve 



dS x 



dx 

S X &Z OD z 



(40) 



OS* 

dx ' (1 - SI) dx 

b x b z Ob z 



dx (1 - S 2 ) fe 
where the last line follows from 



byb z Ob z 



5,-0,-62) - -^J^-—^- 



d_ 
'dx 



(9i-e 2 ) = 



as. 



OfT-Ds- C7*-5z 



5a; (1 - Si) dx 



(41a) 
(41b) 



noting that 

.4- . 1 1 - f-ir ! 



fc=l,2 



- (-i) 3 - fe s 2 



(42a) 



(1 - s 2 ) 



E(-D 3 - fe i 



fe=l,2 



yl +(-l) 3 -*g, 

- (-i) 3 - fe s 2 

-(zo^i> 2 + b x ), 



1/2 



i(e 3 - fc -e fc ) 



(42b) 



(1 - s 2 ) 

is equal to 

-(1-S 2 ) 1 / 2 J2 e*^ 8 -*-®*): 



fe=l,2 

(-i)V /as, 

2(l + (-l) fe S 2 ) 2 V&r 
1 



/ - + (-l) 3 - kdSz 



dy 



(1-S2) 



(zi>a;0 2 — by)— h (ibybz + Jx)~7; 



(43) 



We combine Eq.'s (|40j) and (|43j to get 

1 



s_L-ve+- 



5r 



SS 2 55 ; 



dS^ + dSy 

dx dy 

(44) 



(1 - S 2 ) V ^ 9y y 9a; 
Finally, notice that 

- a - sir- e e«*-»-*4 (ig + (-D 3 - fc |) 



fc=l,2 



; .-, 5p . Q dp 

l[bx a^ + by d~y 



(45) 



Thus, Eq. (38 1 becomes 

E so = Jf (s ± ■ ve 

1 / as 2 _ as 2 

The last step is to notice that 



(46) 



d 2 r. 



S-Vx S 



1 



(1 - si) v a^ 



Q as, as 



y.r; 



S 2 Sj_ • -R I , 

(47) 



so that 



0= f pn(s ± -v cS +^S-\7xs) d 2 r (48) 



We move on to the second term of Eq. (39), which, by 



since 
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1 /__ S Z R 



(49) 



and thus Eq. (36) yields Eq. 05 



2. Alternative Forms for the Energy Functional iEu 



Equation ( 35 ) can be decomposed into its constituent 



parts, i.e. we can write E = E^e + ^pb + Ej, where 



2 , P 



Eke = J g (VVP) + § (V^) Z d\ 

E PE = / 2 ( Vcff - ^ X r ) 2 + />«(<Su. ' «eff 



^.Vx5)+^(l-l]V d 2 r, 



Ei= (co + aSz + czS 



2^P .2 
a r. 



(50a) 



(50b) 



(50c) 



In the following we give two alternative expressions for 
Epe\ (i) note that 



( VcS ~nx r y = (v cB + K Sx) - (^ 2 (i - s 2 z ) - n 2 

- 2kS± -v cS -2rtxr- v cS , 
so that we are able to write Epe as 

EpE = I 2 (Vcff + KS± ^ + Y S ' WXS 

n X r 



(51) 



•L V9 H 

+ 2 V(S 2 -l) + r 2 ) rf 2 r- 



(52) 



(ii) Similarly, we note that 
{v cff - n x rf = (v cS - ft x r + nS±) 2 - k 2 (1 - S 2 



(2kS± ■ ft x r — 2KV of f • S±) 



(53) 



so that we are able to write E PE as 



Epe=\\ (v efl -nxr + kS ± ) 2 + ^S-WxS 
+ pn£l(-yS x + xS y ) 

+ f(« 2 (S* 2 -l) + (l-tfV) d 2 r. 

(54) 

Some of these formulations of the energy are related to 
some computations in [12) or the hydrodynamic formu- 
lation in 1391. 
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